#include "Element_1D.h"

Point_1::Point_1()
{
    coord = new double [1];
    index = -1;
};

Point_1::Point_1(double _x)
{
    coord = new double [1];
    coord[0] = _x;
    index = -1;
};

Point_1::Point_1(double _x, int _i)
{
    coord = new double [1];
    coord[0] = _x;
    index = _i;
};


Point_1::~Point_1()
{
    delete [] coord;
};

int Point_1::get_index() const
{
    return index;
};

void Point_1::set_index(int _i)
{
    index = _i;
};        


std::ostream& operator<<(std::ostream &_os, const Point_1 &_p)
{
    _os << _p.index << ": " << _p.coord[0];
    return _os;
};


Point_1 operator+(const Point_1 &_p1, const Point_1 &_p2)
{
    Point_1 re(_p1.coord[0] + _p2.coord[0]);
    return re;
};


Point_1::Point_1(const Point_1 &_p)
{
    coord = new double [1];
    coord[0] = _p.coord[0];
    index = _p.index;
};

const Point_1& Point_1::operator=(const Point_1 &_p)
{
    if (this != &_p)
    {
	coord[0] = _p.coord[0];
	index = _p.index;
    }
    return *this;
};

double Point_1::get_coord() const
{
    return coord[0];
};

void Point_1::set_coord(double _c)
{
    coord[0] = _c;
};

Interval::Interval()
{    
};

Interval::~Interval()
{};

Interval::Interval(double _a, double _b)
{
    a.set_coord(_a);
    b.set_coord(_b); 
};

Interval::Interval(const Point_1 &_a, const Point_1 &_b)
{
    a = _a;
    b = _b;
};

void Interval::set_a(const Point_1 &_p)
{
    a = _p;
};

void Interval::set_b(const Point_1 &_p)
{
    b = _p;
};

const Point_1& Interval::get_a()
{
    return a;
};

const Point_1& Interval::get_b()
{
    return b;
};

std::ostream& operator<<(std::ostream &_os, const Interval &_i)
{
    _os << "[" << _i.a << ", " << _i.b << "]";
    return _os;
};

Mesh_1::Mesh_1(const Point_1 *_a, int _n_nodes)
{
    interval = new Interval(_a[0], _a[_n_nodes - 1]);
    nodes = new Point_1 [_n_nodes];
    
    for (int i = 0; i < _n_nodes; i++)
    {
	nodes[i] = _a[i];
	nodes[i].set_index(i);
    }
    Interval elements[_n_nodes - 1];
    for (int i = 0; i < _n_nodes - 1; i++)
      {
	elements[i].set_a(_a[i]);
	elements[i].set_b(_a[i + 1]);
      }
    n_grids = _n_nodes - 1;
};

Mesh_1::~Mesh_1()
{
    delete interval;
    delete [] nodes;
};

const Interval& Mesh_1::get_interval()
{
    return *interval;
};


int Mesh_1::get_n_grids()
{
    return n_grids;
};

int Mesh_1::get_n_nodes()
{
    return n_grids + 1;
};

const Interval& Mesh_1::get_element(int _i)
{
    return *elements[_i];
};

Element_1_1::Element_1_1()
{};

Element_1_1::Element_1_1(const Point_1& _a, const Point_1& _b)
{
    a = _a;
    b = _b;
};

Element_1_1::~Element_1_1()
{};

void Element_1_1::set_a(const Point_1& _a)
{
    a = _a;
};

void Element_1_1::set_b(const Point_1& _b)
{
    b = _b;
};

double Element_1_1::phi_a(double _x)
{
    double x0 = a.get_coord();
    double x1 = b.get_coord();
    return (_x - x1) / (x0 - x1);
};

double Element_1_1::phi_b(double _x)
{
    double x0 = a.get_coord();
    double x1 = b.get_coord();
    return (_x - x0) / (x1 - x0);
};


Quad_Info::Quad_Info()
{
    n_q_points = -1;
};

const double* Quad_Info::get_pAw(int _idx) const
{
    return q_info[_idx];
};


void Quad_Info::set(int _acc, int _n_q_pnts, double **_q_info)
{
    accuracy = _acc;
    n_q_points = _n_q_pnts;
    q_info = new double* [n_q_points];
    for (int i = 0; i < n_q_points; i++)
    {
	q_info[i] = new double [2];
	q_info[i][0] = _q_info[i][0];
	q_info[i][1] = _q_info[i][1];
    }
};

Quad_Info::~Quad_Info()
{
    for (int i = 0; i < n_q_points; i++)
	delete [] q_info[i];
    if (n_q_points >= 0)
	delete [] q_info;
};

int Quad_Info::get_accuracy() const
{
    return accuracy;
};

int Quad_Info::get_n_q_points() const
{
    return n_q_points;
};

Quad::Quad()
{
    quad_info = NULL;
};

void Quad::load_data(const char _file[])
{
    if (quad_info != NULL)
    {
	std::cerr << "Quad table is not empty." << std::endl;
	return;
    }
    std::fstream data(_file);
    data >> n_table;
    quad_info = new Quad_Info [n_table];
    int k = 0;
    while (!data.eof())
    {
	int n_pnts, acc;
	double **q_info;
	data >> acc >> n_pnts;
	q_info = new double* [n_pnts];
	for (int i = 0; i < n_pnts; i++)
	{
	    q_info[i] = new double [2];
	    data >> q_info[i][0] >> q_info[i][1];
	}
	quad_info[k++].set(acc, n_pnts, q_info);
	for (int i = 0; i < n_pnts; i++)
	    delete [] q_info[i];
	delete [] q_info;
    }
    data.close();
};

void Quad::list() const
{
    std::cout << "no. of tables: " << n_table << std::endl;
    for (int i = 0; i < n_table; i++)
    {
	int n_pnts = quad_info[i].get_n_q_points();
	std::cout << "table " << i << ":" << std::endl;
	std::cout << "algbra accuracy: " << quad_info[i].get_accuracy() << ", "
		  << "no. of points: " << n_pnts << "." << std::endl;
	for (int ii = 0; ii < n_pnts; ii++)
	{
	    const double* pAw = quad_info[i].get_pAw(ii);
	    std::cout << ii << ", point " << pAw[0]
		      << ", weight " << pAw[1] << std::endl;
	}
    }
};

Quad::~Quad()
{
    delete [] quad_info;
};

int Quad::get_n_table() const
{
    return n_table;
};

const Quad_Info& Quad::get_quad_info(int _acc) const
{
    int i = 0;
    bool found = false;
    for (; i < n_table; i++)
	if (quad_info[i].get_accuracy() > _acc)
	{
	    found = true;
	    break;
	}
    if (!found)
	std::cerr << "The accuracy is not satisfied." << std::endl;
    return quad_info[i]; 
};

double Quad::l2g(double _a, double _b, double _t) const
{
    return (_t + 1.0) * (_b - _a) * 0.5;
};

double Quad::g2l(double _a, double _b, double _x) const
{
    return (-1.0 + 2.0 / (_b - _a) * _x);
};

double Quad::Jacobian(double _a, double _b) const
{
    return (_b - _a) * 0.5;
};

